home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / graphic / frasr182.zip / LYAPUNOV.ASM < prev    next >
Assembly Source File  |  1993-04-07  |  20KB  |  495 lines
  1. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  2. ; Wesley Loewer's attempt at writing his own
  3. ; 80x87 assembler implementation of Lyapunov fractals
  4. ; based on lyapunov_cycles_in_c() in miscfrac.c
  5. ;
  6. ; Nicholas Wilt, April 1992, originally wrote an 80x87 assembler
  7. ; implementation of Lyapunov fractals, but I couldn't get his
  8. ; lyapunov_cycles() to work right with fractint.
  9. ; So I'm starting mostly from scratch with credits to Nicholas Wilt's
  10. ; code marked with NW.
  11.  
  12. .8086
  13. .8087
  14. .MODEL medium,c
  15.  
  16. ; Increase the overflowcheck value at your own risk.
  17. ; Bigger value -> check less often -> faster run time, increased risk of overflow.
  18. ; I've had failures with 6 when using "set no87" as emulation can be less forgiving.
  19. overflowcheck EQU 5
  20.  
  21. EXTRN   Population:QWORD,Rate:QWORD
  22. EXTRN   colors:WORD, maxit:WORD, lyaLength:WORD
  23. EXTRN   lyaRxy:WORD, LogFlag:WORD
  24. EXTRN   fpu:WORD
  25. EXTRN   filter_cycles:WORD
  26.  
  27. .DATA
  28. ALIGN 2
  29. BigNum      DD      100000.0
  30. half        DD      0.5                 ; for rounding to integers
  31. e_10        DQ      22026.4657948       ; e^10
  32. e_neg10     DQ      4.53999297625e-5    ; e^(-10)
  33.  
  34. .CODE
  35.  
  36. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  37. ; BifurcLambda - not used in lyapunov.asm, but used in miscfrac.c
  38. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  39. PUBLIC  BifurcLambda
  40. BifurcLambda    PROC
  41. ;   Population = Rate * Population * (1 - Population);
  42. ;   return (fabs(Population) > BIG);
  43.         push    bp              ; Establish stack frame
  44.         mov     bp,sp           ;
  45.         sub     sp,2            ; Local for 80x87 flags
  46.         fld     Population      ; Get population
  47.         fld     st              ; Population Population
  48.         fld1                    ; 1 Population Population
  49.         fsubr                   ; 1 - Population Population
  50.         fmul                    ; Population * (1 - Population)
  51.         fmul    Rate            ; Rate * Population * (1 - Population)
  52.         fst     Population      ; Store in Population
  53.         fabs                    ; Take absolute value
  54.         fcomp   BigNum          ; Compare to 100000.0
  55.         fstsw   [bp-2]          ; Return 1 if greater than 100000.0
  56.         mov     ax,[bp-2]       ;
  57.         sahf                    ;
  58.         ja      Return1         ;
  59.         xor     ax,ax           ;
  60.         jmp     short Done      ;
  61. Return1:
  62.         mov     ax,1            ;
  63. Done:   inc     sp              ; Deallocate locals
  64.         inc     sp              ;
  65.         pop     bp              ; Restore stack frame
  66.         ret                     ; Return
  67. BifurcLambda    ENDP
  68.  
  69. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  70. ; OneMinusExpMacro - calculates e^x
  71. ; parameter : x in st(0)
  72. ; returns   : e^x in st(0)
  73. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  74. OneMinusExpMacro   MACRO
  75. LOCAL not_387, positive_x, was_positive_x
  76. LOCAL long_way_386, long_way_286
  77. LOCAL positive_x_long, was_positive_x_long, done
  78.  
  79. ; To take e^x, one can use 2^(x*log_2(e)), however on on 8087/287
  80. ; the valid domain for x is 0 <= x <= 0.5  while the 387/487 allows
  81. ; -1.0 <= x <= +1.0.
  82.  
  83. ; I wrote the following 387+ specific code to take advantage of the
  84. ; improved domain in the f2xm1 command, but the 287- code seems to work
  85. ; about as fast.
  86.  
  87.         cmp     fpu,387                 ; is fpu >= 387 ?
  88.         jl      not_387
  89.  
  90. ;.386   ; a 387 or better is present so we might as well use .386/7 here
  91. ;.387   ; so "fstsw ax" can be used.  Note that the same can be accomplished 
  92. .286    ; with .286/7 which is recognized by more more assemblers.
  93. .287    ; (ie: MS QuickAssembler doesn't recognize .386/7)
  94.  
  95. ; |x*log_2(e)| must be <= 1.0 in order to work with 386 or greater.
  96. ; It usually is, so it's worth taking these extra steps to check.
  97.         fldl2e                          ; log_2(e) x
  98.         fmul                            ; x*log_2(e)
  99. ;begining of short_way code
  100.         fld1                            ; 1 x*log_2(e)
  101.         fld     st(1)                   ; x*log_2(e) 1 x*log_2(e)
  102.         fabs                            ; |x*log_2(e)| 1 x*log_2(e)
  103.         fcompp                          ; x*log_2(e)
  104.         fstsw   ax
  105.         sahf
  106.         ja      long_way_386            ; do it the long way
  107.         f2xm1                           ; e^x-1=2^(x*log_2(e))-1
  108.         fchs                            ; 1-e^x which is what we wanted anyway
  109.     jmp     done                    ; done here.
  110.  
  111. long_way_386:
  112. ; mostly taken from NW's code
  113.         fld     st                      ; x x
  114.         frndint                         ; i x
  115.         fsub    st(1),st                ; i x-i
  116.         fxch                            ; x-i i
  117.         f2xm1                           ; e^(x-i)-1 i
  118.         fld1                            ; 1 e^(x-i)-1 i
  119.         fadd                            ; e^(x-i) i
  120.         fscale                          ; e^x i
  121.         fstp    st(1)                   ; e^x
  122.         fld1                            ; 1 e^x
  123.         fsubr                           ; 1-e^x
  124.         jmp     done
  125.  
  126. not_387:
  127. .8086
  128. .8087
  129. ; |x*log_2(e)| must be <= 0.5 in order to work with 286 or less.
  130. ; It usually is, so it's worth taking these extra steps to check.
  131.         fldl2e                          ; log_2(e) x
  132.         fmul                            ; x*log_2(e)
  133.  
  134. ;begining of short_way code
  135.         fld     st                      ; x*log_2(e) x*log_2(e)
  136.         fabs                            ; |x*log_2(e)| x*log_2(e)
  137.         fcomp   half                    ; x*log_2(e)
  138.         fstsw   status_word
  139.         fwait
  140.         mov     ax,status_word
  141.         sahf
  142.         ja      long_way_286            ; do it the long way
  143.  
  144. ; 286 or less requires x>=0, if x is negative, use e^(-x) = 1/e^x
  145.         ftst                            ; x
  146.         fstsw   status_word
  147.         fwait
  148.         mov     ax,status_word
  149.         sahf
  150.         jae     positive_x
  151.         fabs                            ; -x
  152.         f2xm1                           ; e^(-x)-1=2^(-x*log_2(e))-1
  153.         fld1                            ; 1 e^(-x)-1
  154.         fadd    st,st(1)                ; e^(-x) e^(-x)-1
  155.         fdiv                            ; 1-e^x=(e^(-x)-1)/e^(-x)
  156.         jmp     SHORT done              ; done here.
  157.  
  158. positive_x:
  159.         f2xm1                           ; e^x-1=2^(x*log_2(e))-1
  160.         fchs                            ; 1-e^x which is what we wanted anyway
  161.         jmp     SHORT done              ; done here.
  162.  
  163. long_way_286:
  164. ; mostly taken from NW's code
  165.         fld     st                      ; x x
  166.         frndint                         ; i x
  167.         fsub    st(1),st                ; i x-i
  168.         fxch                            ; x-i i
  169. ; x-i could be pos or neg
  170. ; 286 or less requires x>=0, if negative, use e^(-x) = 1/e^x
  171.         ftst
  172.         fstsw   status_word
  173.         fwait
  174.         mov     ax,status_word
  175.         sahf
  176.         jae     positive_x_long_way
  177.         fabs                            ; i-x
  178.         f2xm1                           ; e^(i-x)-1 i
  179.         fld1                            ; 1 e^(i-x)-1 i
  180.         fadd    st(1),st                ; 1 e^(i-x) i
  181.         fdivr                           ; e^(x-i) i
  182.         fscale                          ; e^x i
  183.         fstp    st(1)                   ; e^x
  184.         fld1                            ; 1 e^x
  185.         fsubr                           ; 1-e^x
  186.         jmp     SHORT done              ; done here.
  187.  
  188. positive_x_long_way:                    ; x-i
  189.         f2xm1                           ; e^(x-i)-1 i
  190.         fld1                            ; 1 e^(x-i)-1 i
  191.         fadd                            ; e^(x-i) i
  192.         fscale                          ; e^x i
  193.         fstp    st(1)                   ; e^x
  194.         fld1                            ; 1 e^x
  195.         fsubr                           ; 1-e^x
  196. done:
  197. ENDM    ; OneMinusExpMacro
  198.  
  199. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  200. ; modeled from lyapunov_cycles_in_c() in miscfrac.c
  201. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  202. ;int lyapunov_cycles_